Determination of Properties of Composite Materials from the Lamb Wave Propagation: Probabilistic, Interval, and Fuzzy Approaches

One of the most efficient non-destructive techniques for finding hidden faults in a plate is the use of ultrasonic Lamb waves. These waves are reflected by faults, and from this reflection, we can locate the faults. For that, we need to know how the Lamb waves propagate. Their propagation is determined by the dynamic elastic constants C 0 pq , so we must find these constants. These constants cannot be measured directly; instead, we measure the dependence of the speed on frequency c(f), and we must reconstruct C 0 pq from the measured values of c(f). In this paper, we show how this can be done in the presence of probabilistic, interval, and fuzzy uncertainty. 1. Formulation of the Problem Composite materials are very important. Engineers often use composites, i.e., materials obtained by combining several layers of different materials. Composites are becoming main materials for aerospace and submarine applications, where, e.g., an internal layer provides structural strength while the outside layer guarantees protection against the surrounding media. To apply non-destructive testing to materials, we must know their elastic characteristics. Since composite materials are often intended for use in dangerous environments it is important to periodically test the structural integrity of the corresponding structures. For example, we send an ultrasonic wave into the corresponding structure and see whether any faults are causing the wave to differ from its intended path. To reconstruct the faults from the resulting signal, we must know how the sound propagates in this plate. This propagation is determined by the dynamic elastic constants C 0 pq . paragraphWe must determine elastic characteristics from the dispersion curve. In view of the intended use, a reasonable way to measure these constants is by measuring how the ultrasound waves propagate in this composite material, and to reconstruct the values C 0 pq from the results of these measurements. In the absence of faults, a wave travels in a straight line, so the only information we can get is the wave’s velocity c at different frequencies f ; the dependence c(f) is called a dispersion curve. There is a known algorithm c(f; C 0 pq) that describes how c(f) depends on C 0 pq ; see below. We must therefore reconstruct C 0 pq by solving the system of equations c(f; C 0 pq) = c(f) for measured values c(f). There exist methods for reconstructing the elastic constants from the dispersion curve c(f) (see, e.g., [3]), but these methods assume that we can measure the velocity c for frequencies up to 8 MHz. Experiments with such highfrequency ultrasound are very difficult and expensive, so it is desirable to find out what we can compute based on a more easily accessible measurements for f 2 MHz. Why probabilistic, interval, and fuzzy uncertainty. Uncertainty comes from the fact that the values c(f) come from measurements, and measurements are never 100% accurate: there is a random measurement error component, which corresponds to a probabilistic uncertainty; see, e.g., [11]; there are known bounds on a systematic error component, which correspond to interval uncertainty; see, e.g., [2, 4, 6]; in addition to these bounds in which experts are absolutely sure, experts have smaller bounds with reasonable but not absolute certainty; these bounds are naturally described by fuzzy techniques; see, e.g., [1, 5, 10]. 2. Case Study: Monoclinic Transverse Isotropic Material Let us first recall how the dispersion curve depends on the elastic constants; this dependence is described, e.g., in [8] (see also [7]). In elastic materials, usually, stress 0 ij is linearly related to strain e0kl, i.e., 0 ij = 3 X k=1 3 Xl=1 c0ijkl u0kl; (1) where c0ijkl are constants called elastic material constants. From the purely mathematical viewpoint, we can have 3 3 3 3 3= 81 different constants c0ijkl. Due to physical arguments like energy conservation law, we must have a relation between these coefficients, e.g., c0ijkl = c0klij = c0jikl = c0ijlk . A good way to take these relations into consideration is to consider pairs of coordinates (i; j) instead of single coordinates. In general, there are 3 3 = 9 possible pairs, but since c0 does not change if we swap i and j, it is sufficient to only consider pairs for which i j. These are six such pairs: 11, 22, 33, 23, 13, and 12. These pairs are usually numbered so that 11 $ 1, 22 $ 2, 33 $ 3, 23 $ 4, 13 $ 5, 12 $ 6. In these terms, each value c0ijkl is described as C 0 pq , where p corresponds to ij and q corresponds to kl. For example, c01312 becomes C 0 56. Thus, the elastic material constants form a 6 6 matrix. Since c0ijkl = c0klij , we conclude that C 0 pq = C 0 qp, i.e., the matrix C 0 pq is symmetric. In general, we need 6 (6 + 1)=2 = 21 parameters to describe a general symmetric 6 6 matrix. Most materials, however, have additional symmetry that enables us to further restrict the number of parameters. In this paper, we consider the case when the material is monoclinic and transverse isotropic. Monoclinic means that 0x01x02 is a plane of mirror symmetry, and transverse isotropic means that the x01 axis is orthogonal to the plane of isotropy. In this case, we have only five independent constants: C 0 11, C 0 12, C 0 22, C 0 23, and C 0 55. Once we know these constants, we can reconstruct the entire matrix C 0 pq as follows: 0BBBBBBB@C 0 11 C 0 12 C 0 12 0 0 0 C 0 12 C 0 22 C 0 23 0 0 0 C 0 12 C 0 23 C 0 22 0 0 0 0 0 0 C 0 22 C 0 23 2 0 0 0 0 0 0 C 0 55 0 0 0 0 0 0 C 0 55 1CCCCCCCA : (2) Let ' denote the angle between the direction of the wave and the x01 axis. It turns out that it is convenient to use an auxiliary orthogonal coordinate system x1; x2; x3 in which x3 is the same but x1 and x2 are obtained from x01 and x02 by a rotation by the angle '. The corresponding coordinate transformation can be described as xi = 3 Xj=1 ij x0j , where ij = 0@ cos(') sin(') 0 sin(') cos(') 0 0 0 11A : (3) The values of the elastic coefficients in the new coordinates can be determined as cijkl = 3 X i0=1 3 X j0=1 3 X k0=1 3 X l0=1 ii0 jj0 kk0 ll0 ci0j0k0l0 : (4) Let us now describe the formulas for finding c(f). Let d denote the thickness of the plate, the density. Then, for a symmetric wave, the ratio def = ! d 2 c ; (5) where ! def = 2 f is an angular frequency, satisfies the equation D11 G1 cot( 1) +D13 G3 cot( 3)+ D15 G5 cot( 5) = 0; (6) where G1 def = D23 D35 D33 D25; (7) G3 def = D31 D25 D21 D35; (8) G5 def = D21 D33 D31 D23; (9) D1q def = C13 + C36 Vq + C33 q Wq ; (10) D2q def = C55 ( q +Wq) + C45 q Vq ; (11) D3q def = C45 ( q +Wq) + C44 q Vq ; (12) Vq def = K11( q) K23( q) K12( q) K13( q) K13( q) K22( q) K12( q) K23( q) ; (13) Wq def = K11( q) K23( q) K12( q) K13( q) K12( q) K33( q) K23( q) K13( q) ; (14) K11( ) def = C11 c2 + C55 2; (15) K12( ) def = C16 + C45 2; (16) K13( ) def = (C13 + C55) ; (17) K22( ) def = C66 c2 + C44 2; (18) K23( ) def = (C36 + C45) ; (19) K33( ) def = C55 c2 + C33 2; (20) and q are the roots of the equation 6 +A1 4 +A2 2 +A3 = 0; (21) ordered in such a way that 2 = 1; 4 = 3; 6 = 5; (22) and A1 def = C11 C33 C44 C2 13 C44 +2C13 C36 C45 2C13 C44 C55 + 2C13 C2 45 2C16 C33 C45+ C33 C55 C66 C2 36 C55 (C33 C44 + C33 C55 + C44 C55 C2 45) c2; (23) A2 def = C11 C33 C66 C11 C2 36 2C11 C36 C45 C11 C44 C55 C11 C2 45 C2 13 C66+ 2C13 C16 C36 + 2C13 C55 C66 C2 16 C33 + 2C16 C36 C55 (C11 C33 + C11 C44 C2 13 2C13 C55) c2+ ( C16 C45 + C33 C66 C2 36 2C36 C45) c2+ (C44 C55 C2 45 + C55 C66) c2+ (C33 + C44 + C55) 2 c4; (24) A3 def = C11 C55 C66 C2 16 C55 (C11 C55 + C11 C66 + C2 16 + C55 C66) c2+ (C11 + C55 + C66) 2 c4 3 c6; (25) def = C33 C44 C55 C33 C2 45: (26) 3. Forward Problem: Auxiliary Algorithm Why forward problem. Our objective is to find the elastic constants C 0 pq that explain the observed dispersion curve c(f). To achieve this objective, we must be able, given the constants C 0 pq , to determine the corresponding dispersion curve. In terms of applied mathematics, in order to solve the inverse problem of finding C 0 pq from c(f), we must first be able to solve the forward problem: given C 0 pq , find c(f). Even solving this forward problem is not easy because the above formulas only provide us with an implicit relation between C 0 pq and c(f). It is, however, possible to extract an algorithm from the above formulas. Natural idea. A seemingly natural idea is, given f , to determine c for this f . Then, when we repeat the corresponding computations for different values f = f 0, f = f + f , f = f + 2 f , . . . , f+ = 2 MHz, we get the desired dispersion curve. This was the first algorithm we implemented. Limitations of the natural idea. It turns out that, due to implicitness of the above relations, this algorithm requires solving a lot of implicit equations and runs for hours on a PC. New idea. It turns out that we get a much faster algorithm for reconstructing a dispersion curve (which runs for minutes instead of hours) if we start with different values of velocity c = c , c = c + c, . . . , c = c+, and look for frequencies that correspond to these different values of velocity. For f 2 MHZ, the velocity range [c ; c+] is usually between 1 and 4 km/sec. Input to the forward algorithm. The input to the new algorithm consists of: the constants C 0 11, C 0 12, C 0 22, C 0 23, and C 0 55; the density ; the angle '; and the velocity c. New algorithm: description. The algorithm consists of the following steps: 1. compute ij by using the formula (3); 2. compute all the elements of the matrix C 0 pq by using the formula (2); 3. compute the values Cpq by using the formula (4); 4. compute by using the formula (26); 5. compute the values A1, A2, and A3 by using the formulas (23)–(25); 6. solve the equation (21); we use NSolve program from the Mathematica package; this program produces six roots r1; : : : ; r6; 7. order the roots so that they satisfy the property (22); for this, we do the following: as 1, we take r1; as 2, we take r1; as 3, we take the first of the values ri that is different from r1 and r1; as 4, we take 3; as 5, we take the first of the values ri that is different from 1 and 3; as 6, we take 5. 8. for each of the six roots 1; : : : 6, we compute the following: we compute Kmn( q) by using the formulas (15)–(20); we compute Vq and Wq by using the formulas (13)–(14); we compute Diq by using the formulas (10)– (12); 9. then, we compute Gi by using the formulas (7)–(9); 10. next, we solve the equation (6) with the unknown ; for this, we also use a Mathematica equation solver; 11. based on , we reconstruct ! as (2c )=d and then f as !=(2 ): Interpolation leads to a dispersion curve. Once we know the values f(c ), f(c + c), . . . , f(c + k c), . . . , corresponding to different c, then, for a given f , we can find c(f) as follows: first, by using bisection, we find k for which f(c + k f) f f(c + (k + 1) c); (27) we can now use the linear interpolation on the interval [f(c + k f); f(c + (k+1) c)] and determine the value c corresponding to f as follows: – for f = f(c + k f), we have c(f) = c + k f); – for f = f(c + (k + 1) f), we have c(f) = c + (k + 1) f); therefore, c = c + k c+ c f f(c + k c) f(c + (k + 1) c) f(c + k c) : (28) 4. Inverse Problem Under Probabilistic Uncertainty: Formulation and Challenges Formulation of the problem in precise terms. We measure the velocity for different values of frequency f1; f2; : : :; based on the resulting values ck = c(fk), we must reconstruct the parameters C 0 pq . In the probabilistic setting (see, e.g., [11]), a natural idea is to use the least squares method, i.e., to look for the values C 0 pq for which S def = Xk (ck c(fk; C 0 pq))2 ! min C0 pq ; (29) where c(f; C 0 pq) denotes the result of the forward algorithm. In addition to these “most probable” values C 0 pq , we may also want to describe all the values that are consistent with the given measurement results, i.e., all the values C 0 pq for which S =Xk (ck c(fk; C 0 pq))2 ; (30) where the constant depends on the measurement accuracy, on the number of measurements, and on the desired confidence level. This problem is difficult to solve. There are known methods for solving optimization problems, starting with the gradient techniques. It turns out that due to a complex character of the dependence c(f; C 0 pq), the function S has a lot of local minima, so gradient method gets stuck in them, without going to the desired global minimum. A known modification of gradient descent, specifically designed for such situations, is to start several times with a random starting point, and find the smallest of the resulting (local) minima. Alas, the dependence c(f; C 0 pq) is so complex that this modification requires up to 10,000 iterations– and takes several days to run. Yet another alternative is to apply exhaustive search: choose grid so that the change from a point to the next is smaller than the measurement error (this way, we do not miss the correct value), and test all values C 0 pq from the grid. The main problem with this approach is that is still takes too much time. 5. Inverse Problem Under Interval and Fuzzy Uncertainty: Formulation Case of interval uncertainty. For interval uncertainty (see, e.g., [2, 4, 6]), the only information that we have about the error with which we measure velocity c is the upper bound on this measurement error. In this case, it is natural to look for the values C 0 pq for which, for every measurement k, jck c(fk; C 0 pq)j : (31) Comment on interval uncertainty. In the probabilistic approach, we know the probability of different possible values of the measurement error. As a result, it makes sense to select a single set of values C 0 pq as the “most probable” ones. In case of interval uncertainty, we do not have any information about the relative frequency of different values of measurement error. So, all we can do is return the interval of possible values of each of the elastic constants, without specifying any values within these intervals. Case of fuzzy uncertainty. Fuzzy uncertainty corresponds to the case when, in addition to the upper bound on the measurement error, experts provide us with several bounds = 1 > 2 > : : : > m that correspond to decreasing levels of their certainty: an expert is 100% certain that the measurement error is bounded by ; with a certain degree of certainty 2 < 100%, an expert is sure that the error cannot exceed 2; with a somewhat smaller degree of certainty 3 < 2, an expert is sure that the error cannot exceed 3; etc. In this case, for each possible set of elastic constants C 0 pq , we can determine the degree = 1 with which these values are possible as = i, where i is the largest integer for which jck c(fk; C 0 pq)j i (32) for all measurement results ck. Thus, instead of simply providing a user with the set of all possible combinationsC 0 pq , we provide the user, for each combination, with a degree with which this combination is possible. In other words, instead of a (crisp) set of possible combinations, we get a fuzzy set. To be more precise, we get a nested family of sets corresponding to different levels of certainty; such a nested family of sets is indeed equivalent to a more traditional definition of fuzzy set [1, 5, 10]: if a traditional fuzzy set is given, then different sets from the nested family can be viewed as -cuts corresponding to different levels of uncertainty . 6. Solution: Probabilistic and Interval Cases What can we reconstruct? An experimental analysis. Before designing algorithms for solving these problems, we first decided to check whether it is possible to reconstruct all the elastic constants C 0 pq from the measurements of the dispersion curve c(f) limited to f 2 MHz. To check that, for each angle ', we fixed the values of all the constants but one, changed the selected C 0 pq from its smallest possible valueC pq to its largest possible valueC+ pq , and checked how much this change affects c(f). It turns out that for each angle, for some of the constants, the resulting change in c(f) is less than 2%. Since the most accurate measurements provide an accuracy of 2%, this means that by measuring c(f) for f 2 MHz, we cannot reconstruct the corresponding constants C 0 pq . As a result, for each angle, instead of trying to find all five elastic constants, we fix the values of those that we cannot determine, and only try to find the values of the remaining constants. Basis for this solution: sensitivity analysis. In the process of determining whether each elastic constant C 0 pq affects the dispersion curve, we not only decide which constants affect c(f) and which constants do not. For those elastic constants that do affect c(f), we also get some numerical information on their effect, i.e., an information on how sensitive c(f) is to each of these constants. Before solving the inverse problem, we extended this sensitivity analysis as follows. From the literature, we selected several representative combinations of elastic constants. Then, for each angle ', and for each constant C 0 pq , we did the following: for each representative combination r, we computed the values c(f) with this combination and the values c0(f) when the value of the selected constant is increased by 10%; based on these computations, we computed the guaranteed relative effect of this constant as dr;pq def = max k jc0(fk) c(fk)j c(fk) (33) for the interval case and dr;pq def = vuut 1 #k Xk c0(fk) c(fk) c(fk) 2 (34) for the probabilistic case; we took the largest of these values: dpq def = max r dr;pq . The meaning of dpq is simple: if C 0 pq changes by 10%, then c(f) changes by dpq ; similarly, if C 0 pq changes by, say, 5%, then c(f) changes by 0:5 dpq , etc. In particular, for each variable C 0 pq , if we go from the center of the full range [C pq ; C+ pq ] of its possible values to the entire range, we can determine the largest possible relative change cpq in c(f). If we allow all constants to change, then the largest possible change in c(f) cannot thus exceed c def = Xpq cpq . Resulting algorithm. This algorithm is based on the branch and bound idea well known in numerical optimization; see, e.g., [2, 4, 6]. We want to find values within the box-shaped domain [C 11; C+ 11] : : : [C 55; C+ 55]: At each step of this algorithm, we will have a set of boxes that are guaranteed to contain all possible solutions C 0 pq . Initially, we have only one box–the original domain. At each step, we first take each box and divide each side of this box in half (“bisect” each box). As a result, for the case when we can determine all 5 constants, we get 25 = 32 new boxes; for the case when we can only determine 3 variables, we get 23 = 8 new boxes, etc. For each new box, we find its center C 0 pq , and apply the forward algorithm to find the values c(f) corresponding to this center. After k bisections, the size of each resulting box along each axis C 0 pq is exactly 2 k of the original size; thus, the effect of change of C 0 pq within this box cannot exceed the 2 k of the original change, i.e., is 2 k cpq . Hence, when all variables are allowed to change within the box, the resulting change in c(f) cannot exceedP(2 k cpq) = 2 k c. Thus, if the central value of c(f) differs from the observed dispersion curve by > 2 k c, this means that none of the values from this box are possible and thus, this box can be dismissed. We repeat this elimination procedure until the error 2 k c corresponding to the box size becomes smaller than the measurement error (2%); then, all resulting central points C 0 pq are returned as possible values of elastic constants that are consistent with measurement results. On a standard PC, the new method takes hours instead of days. 7. Solution: Fuzzy Case Natural idea. In fuzzy case, instead of a single error bound , we have several error bounds corresponding to different certainty levels i. A natural idea is therefore to apply the above algorithm to each of these levels. As a result, for each level, we get the set of all combinations C 0 pq that are, within this level, consistent with the measurement results. In other words, for each level, we get an -cut of the desired fuzzy set. Limitations of the natural idea. The main problem with this idea is that to handle the fuzzy case, we need to repeat the interval algorithm as many times as there are different fuzzy levels. For a reasonable case of 10 levels, we thus need 10 times more computations time than for the interval case. Thus, if an interval algorithm required hours, we need days. This is too long. How can we speed up these computations? New idea. The possibility to decrease the computation time is based on the fact that the level sets are nested; thus, if we already know that a point belongs to a set, we know for sure that it belongs to the larger level set. Thus, we can speed up computations as follows: In the above interval algorithm, at each iteration, in addition to dismissing the boxes, we mark the remaining ones by the corresponding level , and update these markings as we decrease the box size.